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MAPPING CLOUDS AND TERRAIN OF EARTH-LIKE PLANETS FROM PHOTOMETRIC VARIABILITY: 
DEMONSTRATION WITH PLANETS IN FACE-ON ORBITS 

Hajime Kawahara 1 and Yuka Fujii 2 
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ABSTRACT 

We develop an inversion technique of annual scattered light curves to sketch a two-dimensional 
albedo map of exoplanets in face-on orbits. As a test-bed for future observations of extrasolar ter- 
restrial planets, we apply this mapping technique to simulated light curves of a mock Earth-twin at 
a distance of 10 pc in a face-on circular orbit. A primary feature in recovered albedo maps traces 
the annual mean distribution of clouds. To extract information of other surface types, we attempt to 
reduce the cloud signal by taking difference of two bands. We find that the inversion of reflectivity 
difference between 0.8-0.9 and 0.4-0.5 fim bands roughly recover the continental distribution, except 
for high latitude regions persistently covered with clouds and snow. The inversion of the reflectivity 
difference across the red edge (0.8-0.9 and 0.6-0.7 jum) emphasizes the vegetation features near the 
equator. The planetary obliquity and equinox can be estimated simultaneously with the mapping 
under the presence of clouds. We conclude that the photometric variability of the scattered light will 
be a powerful means for exploring the habitat of a second Earth. 

Subject headings: astrobiology - Earth - scattering - techniques: photometric 



1. INTRODUCTION 

Recent discovery of rocky exoplanets provides hope 
for opening new doors for research of exoplanets whic h 
harbor life (e.g. iLeger et alJ 12009b iBatalha et all 1201 lh . 
Spectral bio markers such as oxyg e n, ozone, water, 
meth ane (e.g. iDes Marais et al.ll2002t iKaltenegger et al.1 
I2010L and references therein) and a ch aracteristic feature 
of plants known as the red edge (e.g. ISeager et al.l 120051 : 
iKiang et al.l l2007bl laj) will play a key role in searching 
for life on habitable planet candidates. It will also be 
important to understand the environment of the plane- 
tary surface, in other words, the habitat of the planet. 
Indeed, diverse surface environments on the Earth in- 
cluding continents, ocean, and meteorological condition 
serve as the backbone of biodiversity. One of the promis- 
ing approaches to know the landscape of the terrestrial 
exoplanets is to identify surface components using the 
scattered light of the planets through the direct imag- 
ing observations. In this field, the Earth itself has been 
considered as a useful test-bed for fu ture investi g ations 
for extrasolar terrestrial exoplanets. iFord et all (j2001[ ) 
demonstrated that the cloud distribution and inhomoge- 
neous surface of the Earth generate photometric variabil- 
ity of scattered lights due to spin rotation. Thereafter, 
several authors atte mpted light curve inversion to charac- 
terize surface types (ICowan et al.l 120091: lOaklev fc Cash! 
I^MlFuvTTet alJl20int iCowan et al.ll2011l) . 

One concern with the diurnal inversion is separation 
of clouds and other surface components since clouds can 
have an enormous contribution to total reflection light (~ 
70 % in energy for the Earth). The diurnal mapping with 
the single-band phot ometry is significantl y polluted by 
the cloud reflection (lOaklev fc" Cash 2009). Diurnal in- 
version techniques with multiband photometry have been 
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proposed to separate surface types, through t he PCA 
analysis with 7 bands ()Cowan et al.ll2009l [201 If) and the 
albedo modeling with 5 bands (jFuiii et al.ll2011f) . 

The surface area facing on the observer and illuminated 
by the host star changes according to spin rotation and 
orbital revolution in different ways. We demonstrated 
that the two-dimensional mapping using the diurnal and 
annual multi-band photometric variations is possible for 
a clou d-free Earth under the a ssumption of the albedo 
model (jKawahara fc Fuiii|[2010l hereafter KF10). For the 
inversion with the diurnal and annual light curves, the 
cloud signal could be serious because season al variation 
of clo uds is larger than diurnal variation (e.g. lPalle et al.l 
2008). In this paper, we improve our previous method 
of the two-dimensional mapping to be applicable to a 
cloudy planet by using a single band or a few bands with- 
out any assumptions of albedo models. This problem is 
basically ill-posed and needs inversion technique used in 
the field of tomography. Hence we name this mapping 
method "spin-orbit tomography". We test this method 
by applying to light curves of a mock cloudy Earth-twin 
based on real satellite data. 

Another virtue of the two-dimensional mapping is the 
planetary obliquity measurement as discussed in KF10 
for the cloud-free case. The obliquity is one of the few ob- 
servables that have potential to constrai n the formation 
scena rio of the terrestrial planets (e.g. K okubo fc Idal 
2007). We also discuss the obliquity measurement for 
a cloudy planet. 



2. SPIN-ORBIT TOMOGRAPHY 

In general, scattered lights from the planetary surface 
are characterized by the bidirectional reflectance distri- 
bution function (BRDF), which is a function of three an- 
gles among the incident ray, the line of sight, and the nor- 
mal direction of surface. Our inversion assumes the Lam- 
bert (isotropic) reflection to approximate real BRDFs of 
the surface. On this approximation, the intensity of the 
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scattered light is expressed as 



h = 



R 2 F b 



W(^0,e,*;C,e s )m(0,0;6) sin 



5vi 



(1) 



where R is the planetary radius, Fb is the incident flux of 
photometric band b, and m(<j>, 9; b) is the albedo of band 
6 at the planetary colatitude 9 and longitude (p. In this 
paper, we only consider a face-on circular orbit since the 
anisotropic scattering does not affect the phase curve. 
Figure Q] displays the schematic configuration of the sys- 
tem. We denote the spin motion by the angle between 
(<j) = 0,9 = tt/2) and the direction of equinox, 3> and or- 
bital revolution by the ecliptic longitude of the host star, 
measured from the direction of equinox. Planetary 
obliquity £ and the ecliptic longitude at the reference 
time to, ®s specify the spin axis. The 9s is related to 
Equinox Day as t oq — to — (27m — 0s)/w or b where w or b is 
the orbital angular velocity. The weight function is de- 
fined as W((j>, 0; ©; C, ©s) = (es-en)(eo-e R ) where e s , 
eo, and en are unit vectors from the planetary surface 
to a host star, from the surface to an observer, and from 
the planetary center to the surface, respectively. The in- 
tegration is performed over the illuminated and visible 
area Syi defined as the region satisfying es • &r > and 
eo • e R > 0. Taking the coordinate as shown in Fig- 
ure [TJ one obtains eo = (0,0,1), es = (cos 8, sin 0, 0) 
and en — (cos (<f> + $) s'm9, cos ( sin (0 + $) sin 9 + 
sin ( cos 9, — sin ( sin (</> + $) sin 9 + cos ( cos 9) . We set 
R = 1 though an observable is R 2 m((f>,9; b). We use the 
light curve normalized by Fb, d(t;b) = nlb(t)/ Fb. In re- 
ality, Fb can be estimated from the orbital distance of 
the planet, the observed host star flux, and the distance 
between the system and the Earth. 




t = t * 

Fig. 1. — Schematic configurations of the planetary system in 
a face-on orbit. The orbit lies in the x-y plane, corresponding to 
the directions of the equinox (x-axis) and the solstice (j/-axis) of 
the planet. The z-axis is defined by both the normal of the orbital 
plane and the direction of observers. The obliquity £ is defined 
as the angle from the z-axis to the spin axis and the line of sight, 
respectively. The ecliptic longitude measured from the equinox 
is denoted by 0. We define ©5 by the ecliptic longitude at the 
reference time t = to so as to describe = u OI ^(t — to) + 0g. 



By discretizing the planetary surface (<j>j,9j), the light 
curve di = d(ti ; b) with observational noise Si (i = 



1, 2, N) can be modeled as 

d^G(C,Q s )m + e, (2) 

rw-(&,^>6(*i).*(*i);C,es)A«« 

[ elsewhere, 

where Aw s and rrij = m{(f>j, Of, b) are the solid angle and 
albedo of the j— th surface pixels (j — 1,2,..., M). In 
this paper, we use t he HEALPix with M = 768 for the 
surface pixelization (|G6rski et al.ll200"5t ). 

We simultaneously solve m, C , and 9<? from d using 
the Tikhonov regula rization (e.g. Mcnkc 1989; iTarantolal 
2005; Hansen 2010), which balances between observa- 
tional noise and spatial resolution of the surface. The 
misfit function of the Tikhonov regularization is given 
by 
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(3) 



where <7j is the standard deviation of the noise. We adopt 
the average of data as a prior of the model rhi = (db). 
From the Bayesian viewpoint, the minimization of Q\ is 
regarded as maximization of a posteriori likelihood as- 
suming a Gaussian prior for m with the average rh and 
the covarian ce matrix A~ 2 J, and an uniform prior for £ 
and @ s (e.g. lTarantolall2005( ). The solution which mini- 
mizes Q\ for given £ and 0s is expressed as 

(4) 
(5) 



m cst , A = VTi\U (d~Grh) + rh 



K 2 + 



The 



where we define c?j = dj/cj and = Gy/(Tj 
orthogonal matrices U and V are given by the singu- 
lar value decomposition of G = UAV T and n t is the 
i-th eigenvalue of G, th at is, the z-th component of the 
diagonal matrix A (e.g. iHansenl I2010D . The 8ij is the 
Kronecker delta. We search for the best-fit £ and ©s 
by the Nelder-Mead method for different A. Finally, the 
regularization parameter A is determined by the L-curve 
criterion, which is the maximum curvature point of the 
mod el norm jm es t a — rh\ versus residuals \d — Gm es t,\\ 
plot (|HanseiJl20iq) . 

3. RESULTS 

Here we demonstrate the spin-orbit tomography by 
applying it to simulated light curves of an Earth-twin. 
The light curves were comput ed with a line-by-line radia- 
tive transfer cod e RSTAR6B ()Nakaiima k Tanakalfl98l 
iNakaiinial Il983f ) in combination with empirical data of 
monthly land reflectivity, monthly snow cover fraction, 
and daily cloud dis tribution in 2008 wh ich are provided 
by Terra/MODIS (jDorothv et al.|[200l . The effect of 
multi-scattering, molecular lines, and the anisotropy of 
the scattering are included by RSTAR6B. We put cloud 
layers made of pure liquid water at an altitude of 4-6 km. 
We note that the different altitude assumptions (0-3 and 
7-20 km) do not change the reflectivity of clouds in the 
three bands significantly (<5 %). The cloud property 
is specified by two parameters, cover fraction and opti- 
cal thickness. The reflection by oc ean is modeled b y the 
Fresnel reflection of a wavy ocean (lNakaiimalll983l) . The 
detailed description is found in lFuiii et al.l ( 20111 ). 
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Fig. 2. — The panel a displays the mock light curve of the blue band (0.4 — 0.5/^m) of an Earth twin at a distance of 10 pc with a 4 
m telescope and a hypothetical detector. The data during 366 days are stacked for 6 days. Hence there are 30 X 61 = 1830 data points. 
Wider versions of two shaded regions are inserted in the mini panels. The predicted curve is drawn by solid lines. The recovered maps on 
the Mollweide projection are shown in the panel b (blue: 0.4 — 0.5^m), c (orange: 0.6 — 0.7^m), and d (NIR: 0.8 — 0.9/im). Magenta curves 
outline the boundary of continents and oceans. We fix £ and ©g to the input values 90° and 0°, respectively. The panel e displays annual 
mean of optical depth of clouds computed from the MODIS level 3 Daily Joint Aerosol/Water Vapor/Cloud Product in 2008 (Dorot hy et alj 
2006). The broad bands of clouds at high latitudes correspond to the polar front and the narrow band of clouds at the equator is known 
as Intertropical Convergence Zone. 
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Fig. 3. — Representable reflection spectra of clouds computed by 
RSTAR6B (gray), soil (magenta), vegetation (green), snow (cyan) 
and water (blue) taken from ASTER spectral library. The Entisol- 
Ustifluvent, Inceptisol-Dystrochrept, and Mollisol-agialboll corre- 
spond to "Brown to dark brown silt loam " , "Dark yellowish brown 
micaceous loam" , and "Dark grayish brown silty loam" , respec- 
tively. Shaded regions indicate three bands we use. 

In this paper we focus on photometric variability in 
three bands: b = NIR (0.8-0.9 /im), orange (0.6-0.7/xm), 
and blue (0.4-0.5^m) since these bands have large differ- 
ences for soil and vegetations. These bands correspond 
to three of the bands of the EPOXI mission, which wa s 
used for validation of the simulation ([Fuiii et al.ll2011l) . 
For demonstration purpose, we assume a face-on circu- 
lar orbit and adopt £ = 90°, which is the best geometry 
for a face-on orbit (KF10). Though this is not the case 
for general inclination, we postpone this problem until a 
forthcoming paper. 

We computed the synthetic scattered light of the Earth 
for 30 epochs per day (in 0.8 hour interval) over 366 days 
and stacked them for 6 days (N — 61 x 30) . Since the spin 



rotationa l period can be m easured through periodogram 
analysis (jPalle et alJl2008l ). we assume that the w sp i n is 
known as well as w or b. The Gaussian noise is imposed so 
as to mimic the observation with signal-to-noise ratio of 
20 for each bin during 0.8 x 6 = 4.8 hours, corresponding 
to a one year observation of an Earth twin at a distance 
of 10 pc with a 4 m telescope and a hypothetical detec- 
tor (KF10). We estimate observational noises assuming 
a future mission with the occulter s ystem, such as the 
occulting ozone observ atory (03; e.g. (Kasdin etail [20101: 
ISavranskv et al.ll2010D . including the exzodiacal light (23 
magarcsec 2 ), dark noise (10 _3 counts s^ 1 ), read noise (2 
Vcounts/read), quantum efficiency (0.91), and end-to- 
end effic iency (0.5) with the p arameters used in KF10 
(see also ISavranskv et"aLll2010t ). While the starlight sup- 
pression ratio is difficult to discuss in general, a t ypical 
supp ression ratio of an oc culter is below 10~ 10 (e .g. ICashl 
[20p), even 10~ 12 for 03 (ISavranskv et al.lf2010h . Other 
noises we assume are comparable to planet's reflection, 
which corresponds to a 10" 10 order of magnitude of the 
starlight. Therefore we simply ignore the star suppres- 
sion noise. Figure [2] a shows the resultant light curves of 
the blue band for £ = 90.0°. 

We perform the inversion of the light curve for each 
band. While we take £ and ©s as fitting parameters, 
the estimated values from the single-band photometry is 
systematically biased a listed in Tabled] (blue). This is 
likely due to the anisotropic scattering and variation of 
clouds. Hence we fix £ and 6s to the input values in 
Figure [2] We will discuss about the obliquity determina- 
tion later. Figure [2] b-d displays the recovered maps of 
the blue (panel b), orange (c), and NIR (d) bands. The 
recovered maps in all bands exhibit a primary feature 
of high (low) reflectivity at high (low) latitude. Com- 
paring with the annual mean of the cloud optical depth 
( Fig. [2] e), one can interpret the primary feature as 
the mean cloud distribution. While the spatial resolu- 
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tion of the inversion is too poor to recover the narrow 
band of persistent clouds seen at the equator, known as 
the Intertropical Convergence Zone, broad bands of the 
polar front cloud at high latitude is well recovered in 
these maps. While short-time scale variations of clouds 
in weeks or months make systematic residuals between 
the prediction and data, as indicated in Figure [5] a, these 
variations do not affect recovered mean features of clouds 
so much. 



20061) and simulations dTinetti et al.l l2006allbl ; 
Montahes-Rodriguez et all l2006t lArnold et al.l T2009). 



TABLE 1 
Estimated obliquity ( and 0g. 



Input 



Estimated 



C = 90°, S = 0° 



blue (S/N=20) 

C = 57.3° ± 1.0° 



12.6° ± 0.9° 



C = 9o°, e s = 

C = 60°, S = 
C = 45°, S = 
C = 30°, 6 S = 
C = 23.4°, e s 



NIR 

0° 
0° 
0° 
0° 
= 0° 



blue (S/N=20) 

= 90.0° ± 1.1°, © s = 

= 59.7° ± 3.1°, S 

= 50.4° ± 2.5°, S 

= 38.1° ±3.0°, S 

= 30.1° ± 2.3°, S 



-1.7° ± 3.2° 
: 15.5° ± 2.9° 
: 17.1° ± 3.9° 
: 19.8° ± 3.3° 
: 23.5° ±3.3° 



NIR - orange (S/N=20) 

C = 90°, 0s = 0° C = 83-5° ± 3.0°, S = -2.0° ± 12.0° 
* 0° < C < 90° and -180° < S < 180° 
The 1 a errors are estimated by the bootstrap resampling. 



Though a primary feature is dominated by clouds, 
slight differences of the recovered maps between different 
bands are due to surface components other than clouds. 
Figure [3] demonstrates several representative examples of 
the surface reflectivity spectra in 0.4 - 1.0 /xm. The re- 
flectivity of soil and vegetation increases as wavelength 
increases. This band dependence causes the differences 
between the recovered maps. Due to almost constant 
reflectivity of clouds (Fig. [3]), the cloud signal is sup- 
pressed by taking the difference between two bands. Ap- 
plying the spin-orbit tomography to the difference vector 
disriR — rfbiue (Fig. [J a), we obtain the NIR-blue map as 
shown in Figure [4] b. The contrast standing out in the 
NIR-blue map traces the approximate continental distri- 
bution of the Earth. The NIR-blue map hardly exhibits 
features at the north America Continent and northern 
part of the Eurasian Continent due to the large cloud 
coverage (Fig. [2]) and almost constant reflectivity of 
snow (Fig. [3]). The reason of negative value on the 
oceans is that Rayleigh scattering dominates the reflec- 
tivity. 

The estimated £ and ©s agree well with the input value 
in this case because the variation and the anisotropic ef- 
fect of clouds are suppressed. We also perform the inver- 
sion for different obliquities listed in Table [TJ The errors 
of £ and ©s are estimated by the bootstrap resampling. 
We find that the inversion of the NIR-Blue bands can 
constrain the obliquity, while the estimate of low obliq- 
uity tends to be difficult due to less available information 
and to be slightly biased. 

We also solve the inversion of difference of two 
bands across the red edge, <Jnir — d, 
[J b). The red edge feature has 
as a potential biomarker by both Earth-shine ob 
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Compared with the NIR-blue map, the NIR-orange 
map (Fig. [4] d) displays larger reflectivity near the 
equator, corresponding to the rain forests in Amazon 
and Southeast Asia. Inhabited exoplanets are likely 
to exhibit the localization of photosynthetic organisms 
because of inhomogeneous insolation and precipitation. 
The enhancement of the edge-like signature at regions 
suited for photosynthesis might support the presence of 
extraterrestrial plants. 

The red edge feature in spatially-unresolved data is 
generally diminished by contributions of oceans and 
clouds. The bar in Figure 0] d indicates the maximum, 
average, and minimum of the reflectivity difference of 
spatially-unresolved light curve (cJnir — <i oran g e ) without 
noise. Since the red edge feature recovered near the trop- 
ical regions is 2-3 times larger than the maximum and 
average of the NIR-orange of the light curve, we con- 
clude that the two-dimensional mapping improves the 
detectability of the red-edge feature. 

In this paper, we have assumed a 4 m tel escope, which 
is sm aller than that of other works (e.g . ICowan et al.1 
2009). For instance, ICowan et al.l (|2009l ) assumed a 16 
m telescope with a coronagraph to achieve 2 % photom- 
etry with 1 hour exposures for an Earth twin at 10 pc. 
Considering the shot noise only, we compute the S/N for 
1 hour exposures and aim telescope: 0.32 and 0.44 
for their and our fiducial (S/N = 5 % for 4.8 hour ex- 
posures) assumptions. Hence, our noise assumption is 
similar to them. The main difference in a telescope size 
requirement comes from that our method uses data in a 
whole year. We also assumed the spin rotati onal period, 
which is a critical parameter for our method. Pall e et al.1 
(2008) suggested that S/N=20 for 0.5 hour exposures m 
a 8-week observation is enough to measure w sp i n by the 
periodogram analysis, corresponding to ~ 3 times larger 
aperture than ours (i.e. ~ 12 m). While a 1 year obser- 
vation might improve statistics, it is unclear whether the 
periodogram works in the same manner when one cannot 
ignore the orbital revolution. Hence the feasibility of the 
spin-orbit tomography in cooperation with the spin ro- 
tation measurement is still under consideration and will 
be discussed in a future work. 

4. SUMMARY 

Using simulated scattered light curves from an Earth- 
twin based on the satellite data, we have demonstrated 
that surface albedo maps can be recovered from annual 
and diurnal photometric variations of terrestrial exoplan- 
ets partially covered with clouds and in face-on orbits. 
We have shown that the inversion of the light curve re- 
covers the approximate cloud distribution. After reduc- 
tion of the cloud signal by taking differences between 
0.8 — 0.9/im and 0.4 — 0.5/xm (or 0.6 — 0.7/xm) , we have 
found that the recovered maps trace continental distribu- 
tion of the Earth. Moreover, the obliquity and equinox 
can be estimated simultaneously with the reflectivity dif- 
ference mapping. 

Though we only consider a face-on orbit in this paper, 
it is essential to extend our technique to general geometry 
for practical use. Especially for low-obliquity planets like 
the Earth, larger area can be recovered in some inclined 
geometry. For the extension to general inclination, the 
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Fig. 4. — The panels a and c display reflectivity differences between the NIR and blue bands and the NIR and orange bands. The 
recovered maps of reflectivity difference are shown in the panel b (NIR-blue) and d (NIR-orange). In this figure we simultaneously estimate 
£ and 0g listed in Table IT1 The black bar in the panel d indicates the maximum (0.02), averaged (0.05) , and minimum (0.08) values of 
the intrinsic reflectivity difference (NIR-orange) of the spatially-unresolved data within a 1 year observation (30 X 366 bins) . 



effect of anisotropic scattering to the phase curve cannot 
be ignored. We will include these effects to the spin-orbit 
tomography and apply it to various geometrical cases 
including transiting planets, low obliquity planets with 
proper inclination, and even tidally locked planets in a 
forthcoming paper (Fujii and Kawahara in preparation). 
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